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0\ ; ABSTRACT 

^^ ■ By direct hydrodynamic simulation, using the Piecewise Parabolic Method (PPM) 

rr . code PROMETHEUS, we study the properties of a convective oxygen burning shell 

1^^ I in a SN 1987A progenitor star ( 20 Mq) prior to collapse. The convection is too 

CN I heterogeneous and dynamic to be well approximated by one-dimensional diffusion-like 
algorithms which have previously been used for this epoch. Qualitatively new 

^ ■ phenomena are seen. 

^^ ■ The simulations are two-dimensional, with good resolution in radius and angle, 

CSJ ■ and used a large (90-degree) slice centered at the equator. The microphysics and 

^I^ . the initial model were carefully treated. Many of the qualitative features of previous 

t~^ i multi-dimensional simulations of convection are seen, including large kinetic and 

^^ I acoustic energy fluxes, which are not accounted for by mixing length theory. Small 

1^ ' but significant amounts of ^■^C are mixed non-uniformly into the oxygen burning 

I ' convection zone, resulting in hot spots of nuclear energy production which are more 

^ • than an order of magnitude more energetic than the oxygen flame itself. Density 

c/3 . perturbations (up to 8 %) occur at the "edges" of the convective zone and are the 

C^ ■ 

• • , result of gravity waves generated by interaction of penetrating flows into the stable 



>$ 



region. Perturbations of temperature and Ye (or neutron excess r]) at the base of 
the convective zone are of sufficient magnitude to create angular inhomogeneities in 
C^ I explosive nucleosynthesis products, and need to be included in quantitative estimates 

of yields. Combined with the plume-like velocity structure arising from convection, the 
perturbations will contribute to the mixing of ^^Ni throughout supernovae envelopes. 
Runs of different resolution, and angular extent, were performed to test the robustness 
of these simulations. 
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1. INTRODUCTION 



Oxygen burning in a convective shell is one of the last, dominant burning stages before core 
collapse of a massive star. Its importance lies not just in the nucleosynthesis occurring by the 
nuclear burning itself, but also in the structure left by convection which might affect subsequent 
supernova evolution ( [Arnett 1969 , Falk fc Arnett 1973| , Arnett et al. 1989). For instance, the 
division between the proto-neutron star and supernova ejecta, known as the "mass cut," is 
expected to appear somewhere near the "edge" of the convective zone. Convection generated 
perturbations in density, velocity, entropy, and Yg (or neutron excess r/) could alter the nature 
of core collapse as compared to spherically symmetric, one-dimensional calculations ( [Arnett 
1977a| , Burrows &i Lattimer 1985| , Bruenn 1985 , Miller, et al., 1991 ; see Arnett 1996| for extensive 



discussion and references). 

The bottom of the oxygen convective shell is expected to be the region where explosive 
nucleosynthesis occurs. Most iron peak isotopes are thought to be produced here. The most 
important individual product of this region would be ^^Ni, which is responsible for powering 
supernova light curves and (after decay to Co) is the source of the 847 and 1208 keV 7-ray lines 
seen in SN1987A (|Leising fc Share 199q , [Tueller et al. 1990| , [Kurfess et al. 1992] , palmer et al. 1993 , 
Ait-Quamar et al. 1992). The explosive O- and Si-burning which produces ^^Ni also produces 



other 7-ray emitting radionuclides such as ^^Ti and ^^Ni (^^Co). Perturbations in Ye, density, and 
temperature, left by convective 0-burning, will greatly affect the relative abundances of individual 
isotopes, due to the sensitivity of nuclear statistical equilibrium (NSE) and quasiequilibrium to 
these parameters ( p^ruran et al. 196^ , Woosley et al. 19"73| , Arnett 1996|) . Studies of NSE have 
shown a neutron excess r] = 1 — 2Ye ~ 0.002 best recreates the relative abundances of iron-peak 
nuclei, with slight changes in 77 altering relative abundances by large factors. A change in YJ, from 
0.499 to 0.4995 changes the relative abundance of ^"^Fe to ^^Fe by roughly a factor of 2. 

We know from observations of SN 1987A that significant mixing of newly synthesized material 
thoughout the envelope must occur ( [Arnett 1988 , Pinto fc Woosley 1988| , Kumagai et al. 1989, 



Graham 1988) . The only physical explanation that agrees with our understanding of massive star 



evolution is that significant density perturbations must appear at the interfaces of composition 
discontinuities and these result in Rayleigh- Taylor (and Richtmeyer-Meshkov) mixing in the ejecta 
( Falk fc Arnett 1973| , Benz &: Thielemann 1990| ). This hypothesis has been tested with various 
fluid hydrodynamics codes by adopting ad hoc distributions of initial perturbation amplitude and 
length scale (|Arnett et al. 1989| , [Fryxell et al. 199l| , pVIiiller et al. 199l| , [Arnett 199i| , [Hachisu et 



al. 1992[ , [Herant fc Benz 1992| , [Herant fc Woosley 199^) . On the whole, calculated ^^Ni velocities 
were in agreement with observations for the bulk of the matter, but not for the highest speeds. 
Not enough Ni could be mixed to velocities as high as 3,000 km s~^, which the observations imply. 
We believe that convective shell 0-burning has a major effect on all of these phenomena. 

In almost all hydrostatic stellar models, convection is treated by enforcing an adiabatic 



gradient, or better, some version of mixing length theory (Prandtl 1925, Vitense 1953, Bohm- 



Vitense 1958; hereafter MLT). The basic tenets of MLT are that (1) convection occurring deep in 
stehar interiors results in an almost completely adiabatic structure, and that (2) excess energy flux 
(above that of radiative diffusion at the adiabatic gradient) is a function of the local superadiabatic 
gradient. The mixing of chemical species is approximated either by complete homogenization, 
or the use of a radial diffusion equation (when MLT formalism is applied). The concept of 
penetrative convection (overshoot from convective regions into adjacent radiative zones) is not a 



natural outcome of MLT, due to the local nature of the theory ( Renzini 1987 ). Various ad hoc 
approaches to overshoot have been applied to stellar evolution; in most the overshoot distance is 
scaled to some multiple of the local pressure scale height. 

Multi-dimensions simulations of turbulent compressible convection have revealed both the 
successes and shortcomings of MLT. One dynamical feature of convection, which has been present 
in every multi-dimensional simulation, is the asymmetry between up and down motions. Upward 
moving elements are typically broader in extent and slower than downward moving elements, 
although the degree of asymmetry is very dependent on the density contrast of the unstable region 
([Graham 19751 Purlburt et al. 1984 |Chan et al. 198^ , [Hurlburt et al. 1986i phan fc Sofia 1986 



Cattaneo et al. 1991). As density contrast increases, so does the asymmetry ([Hurlburt et al, 



1984 ). Pressure fluctuations (the compressible nature of the fluid) can assist in driving descending 



motions, lead to buoyancy braking in ascending motions, and satisfy a Bernoulli relation with 
horizontal velocities ( [Hurlburt et al. 1986| ). The basic notion of buoyancy driven motions 
throughout most of the convective zone has been upheld, but at the boundaries, dissipation and 
compressibilty are the dominant factors in the flow. In the more complex case, in which convective 
regions have nuclear timescales comparable to the mixing timescale, a kinetic (rather than steady 
state) description becomes crucial. 

There are hints that the basic tenet in MLT regarding energy flux still holds true in 
multidimensional simulations. In both 2- and 3-D "large eddy" simulations, in which a sub-grid 
scale turbulence formalism (SGS) has been applied, vertical motions have been seen to breakup 
in less than the simulation dimension (Chan et al. 1982, Chan fc Sofia 19*8^ , Chan fc Sofia 19"89| , 
Hossain fc Mullan 199l| ). In fact, the autocorrelation of vertical velocity along vertical direction is 



a function of the difference in logarithmic pressure, but symmetric with a similar full width-half 



maximum (FWHM) across the entire computational regime ( Chan et al. 1982 , Chan &: Sofia 



1986 ). This means that the convective flux should be a function of the local pressure scale 



height. Amazingly, vertical velocities seem to be correlated over about one pressure scale height 



(Hossain &: Mullan 1991), which is the range used for the past 30 years in stellar evolution models. 
However, even these positive indications are faulty. There seems to be a numerical dependence 
to these results: the correlation is no longer constant over an entire grid when the aspect ratio is 
increased ( phan &: Sofia 1986| ) and the correlation length is dependent on the magnitude of the 
SGS turbulence included in the model. Indeed, in models where SGS turbulent viscosity has been 
omitted in favor of ordinary dynamical viscosity, no breakup of vertical motions has been observed 



([Hurlburt et al. 1984 [Hurlburt et al. 1986|) . 



Penetrative convection has also been modelled in multi-dimensions. In what follows, we 
will use a nomenclature in which penetration refers to mixing (in formally stable zones) that 
is efficient enough to weaken thermal stratification. Otherwise, edge mixing is referred to as 
overshooting ( Massaguer et al. 1984| ). Fast downward plumes are responsible for extensive mixing 
of composition and energy from the unstable region into the lower stable region through gravity 
waves ( Hurlburt et al. 1986| ), although a rapid increase in the stability of adjacent regions reduces 
the ability of plumes to penetrate ( Hurlburt et al. 1994| ). Composition is advected across the 
stability boundary to a "penetration depth," on the timescale for plume transit. In the case of true 
penetration, the region adjacent to the convective zone can be described by two layers: (1) where 
Peclet number > 1, adiabatic structure exists, and motions decelerate due to buoyancy breaking 
and (2) where Peclet number < 1, temperature perturbations damp out from radiative damping 
(Zahn 1991, Hurlburt et al. 1994). Note that the Peclet number is the ratio of the velocity scale v 
times the length scale i, to the thermal conductivity K, Pec = vi/K, and a large Peclet number 



implies ineffective thermal energy transport relative to mass transport; see Landau & Lifshitz 



1959. Our case, in which cooling is by uu emission instead of radiative diffusion, does not seem to 



have been considered by the hydrodynamic community. 

Shell oxygen burning differs from the canonical problem addressed by most previous 



multidimensional convection simulations. The convective regions of ZAMS A and F stars (Latour 
et al., 1981 , Sofia fc Chan 198^ ) and generic convection with ionization ( |Rast et al. 1993 , Rast Sz 



Toomre 1993a , Rast &: Toomre 1993b ) are most comparable. However, such convective regions 
do not contain nuclear energy sources and neutrino energy sinks within the flow. It seems that 
local thermal effects, such as heating from ionization, can significantly affect convective flows. 
The treatment of convection affects neutrino cooled regions, since the efficiency of convection 
determines which adiabat a convective region can achieve. This in turn feeds back upon the 
neutrino emission ( [Aufderheide 1995 , Arnett 1996| ) . 



In previous work (Arnett 1994, Bazan fc Arnett 1994| ), we used our Piecewise Parabolic 
Method fluid dynamics code PROMETHEUS ( [Arnett et al. 1989 ) to begin an examination of the 
problem of the evolution of a 0-burning shell in a 20 Mq, Z = 0.007 star, prior to core collapse. 
Significant density perturbations (of the order of 7%) were found to occur at the boundaries of the 
convective region, along with non-uniform mixing of chemical species within the shell. Here, we 
present more extensive results of models of the same star, with different boundary conditions, grid 
sizes, and grid resolutions, and for a longer evolutionary time. We compare results for different 



boundary conditions to see the effects on the dynamics (Hossain &; Mullan 1993). Computational 
domains are varied in order to see the effects of sound waves on the system. Calculations with 
different grid resolutions are performed to estimate the role of scale lengths on the problem. A 
detailed analysis of the simulations and their astronomical implications is given. 



2. COMPUTATIONAL DETAILS 



2.1. The PROMETHEUS code 

The hydrodynamics code PROMETHEUS is based upon the piecewise-parabohc method 
(PPM, polella &: Woodward 198^ ). The method constructs the physics of the flow between 
grid points by a non-hnear solution of the equations for conservation of mass, momentum, 
and energy (the Riemann problem) rather than the usual Taylor expansion about the grid 
points. This procedure affords greater resolution per grid point, which is highly desirable for 
multi-dimensional problems in stellar astrophysics. Although the computational effort required 
per grid point is higher than other commonly employed numerical methods, the number of grid 
points needed for a given level of accuracy is less (often much less). The net result for PPM 
is a lower total computational effort for given accuracies than competing numerical methods. 
Because the computational load per grid point is greater, a host of physics packages (nuclear 
reactions, radiation, gravity, etc.) may be added before affecting the runtime significantly. Thus 
PPM is well suited for multi-dimensional problems involving significant physics beyond the bare 
hydrodynamics . 

The PROMETHEUS code was originally written by Bruce Fryxell and Ewald Miiller, and 
is an extension of the basic PPM scheme in that (1) it includes an arbitrary number of separate 
fiuids to account for abundances of nuclear species, (2) nuclear reactions and energy generation are 
taken into account, and (3) a realistic (not a gamma law) equation of state is used ( [Colella fc Glaz 



1985). Variations of this original code have been widely employed by other groups (e.g. Blondin & 



Lundqvist 1992, Burrows &: Fryxell 19921 ), including one specifically for the study of convection in 



three dimensions ([Cattaneo et al. 199l| ). PROMETHEUS can support different grid geometries, 
with cartesian, cylindrical, or spherical coordinates in 1, 2, or 3 dimensions as options. There is 
also an option of a moving grid, which allows a 'semi-lagrangian' approach to certain problems 
and adds greater resolution for problems where flows are signiflcant across the entire grid, as in 
explosions. 



2.2. Equation of State 

The equation of state is the sum of components for electrons, ions, and radiation. The electron 
contribution is determined via cubic interpolation in tables for the thermal contributions (including 
e~^e~ pairs), plus degeneracy. Derivatives of the electron pressure and energy as functions of 
temperature and density are obtained from the direct differentiation of the interpolation functions. 
Table entries are logarithmically spaced both in temperature (in degrees Kelvin, Alog(T) = 0.05 
from 7 < log(T) < 10) and in density (in gcm~^, Alog(/?) = 0.125 from 2 < log{p) < 10.5). A 
perfect gas equation of state is smoothly joined for (p,T) regions not covered by the electron EOS 
table, and used for the ion EOS (to which a coulomb correction is added: Landau &: Lifshit^ 



1969). Radiation pressure and energy density are included. 



For the temperatures encountered in quiescent oxygen-burning, nuclear statistical equilibrium 
effects may be ignored, and the effects of quasi-equilibrium seem small. These effects are both 
ignored with regard to the equation of state. 



2.3. Nuclear Reactions 

The algorithms for nuclear reations are similar to those used in prior quasi-static stellar 
evolution studies of massive stars ( [Arnett 1971| ). Twelve species (e~, H, n, ^He, ^^C, ^^O, ^''Ne, 
^^Mg, SiCa, ^^Ni, ^^Co, ^^Fe) are considered for the energy generation and nucleosynthesis 
resulting from pp, CNO, helium (3a), carbon, neon, oxygen, and silicon burning stages. The 
decays of ^^Ni and ^^Co are also included with density dependent decay rates ( Fuller et al. 1985|) . 



The effects of steady state electron captures during oxygen and silicon burning are included; 
Urea processes and their attendant dynamic effects will be dealt with in subsequent work. The 
approximations in the reaction network used to represent oxygen and silicon burning are based 
upon extensive reaction network simulations ( [Thielemann Sz Arnett 1985 , Arnett & Thielemann 
1985i [Arnett 1996|). 



Cooling by neutrino-antineutrino pair emission is included. This process has been identified 
in ID stellar models ( Arnett 1971] ) as competing with nuclear burning in a fashion that determines 



the thermal balance in the oxygen convective region. The thermal conditions of the flame zone, 
in particular, were determined by the bulk of the material in the convective region, which is vV 
cooled, and advected into the flame zone. Using MLT, the heated material would then rise and 
cool by expansion and further vV cooling. Thus, cooling by neutrino-antineutrino pair emission is 
never comparable locally to nuclear heating ( [Arnett 1969| , Arnett 1971 ), rather it is the integral 



effect of heating and cooling across the convective zone which is important for thermal balance. 



We employ the pair emission rates of [Beaudet et al. 1967] , with small revisions for neutral current 
effects ( [Itoh et al. 199^ ). 



2.4. Computational Differences 

There are differences between these simulations and our exploratory work (Arnett 1994) 



Here we map the structure of the star resulting from one-dimensional hydrostatic calculations 
onto the upper and lower boundaries of the grid. Our previous work utilized constant entropy and 
isothermal outflow as the respective lower and upper boundary conditions. We calculate the 2D 
hydrodynamics from the outer region of the Si core to the H envelope, instead of just the O shell 
layer. 



Another difference is the calculation of models with periodic boundary conditions. The initial 
work, as well as some of the simulations discussed here, were calculated with reflective boundary 
conditions, which are known to attract entrained flows and enhance the downward kinetic energy 
flux (Chan & Sofla 1986). Aside from not being realistic for sectors smaller than < < tt, 



reflective boundary conditions can also damp waves of wavelengths larger than the dimension 
of the calculation. Because we do not know exactly how power cascades to smaller scales in 
astrophysical conditions, as compared to box simulations with lower Reynolds' number, we find it 
to be prudent to use periodic boundary conditions rather than reflecting ones. 

We have calculated models over a range of initial grid dimensions, grid resolutions, and 
boundary conditions to see how changes in any of these quantities would change the most critical 
results of this study. Table 1 shows the various descriptive qualities of the models. The differences 
are small (see Appendix). Most of the discussion is illustrated by the high resolution run, case 
D (460x128). A dominant constraint in choice of zoning is that the oxygen flame zone be well 
resolved ( [Arnett 199^ ). 



2.4-1- Numerical Effects 

As a check of the robustness of our calculation, we have also run cases which have the same 
resolution, but a larger or smaller physical dimension. Calculations with a smaller physical domain 
are more sensitive to boundary conditions (and errors), so that larger domains are preferable. 
This is especially true for the oxygen shell burning epoch, because the other parts of the star are 
also evolving. The Si core is treated as an inert lower boundary; numerical experiments (to be 
reported separately) indicate that this approximation is valid over at least the early part of oxygen 
shell burning. The outer boundary is placed near the outer edge of the He core (M^ ~ 6Mq). 
This insures that the acoustic flux is not artificially trapped, includes any back reaction from 
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Table 1: Some Models in this Study. 



^Boundary conditions: R and P denote reflective and periodic boundaries. 



overlying burning shells, and allows us to construct a reasonable two dimensional approximation 
to a SN1987A progenitor model. 

Cases A, B, and D of Table 1 provide a basis in which to examine effects of grid resolution; 
even the coarse case A gave qualitatively similar results. Case E tested the effect of a larger 
angular wedge (135 degrees instead of 90); the results were similar to case B. 

The side boundaries provide a more serious challenge. Reflecting boundaries tend to trap 
waves, so that we have used periodic boundary conditions. Case C is identical to case B except 
for the use of reflecting instead of periodic boundaries; the effect is small for our simulations. 

As an attempt to understand possible geometric peculiarities caused by the choice of grid 
geometry, we have used both an {r,9) coordinate system (Table 1), and some cases which used 
an (r, (/>) coordinate system ( Bazan Sz Arnett 1994| ). The qualitative nature of the results is 
unaffected. 



3. DISCUSSION 

Convection in oxygen shells of massive stars is not driven in quite the same manner as 
convection in stellar envelopes or terrestrial phenomena. The equations to be solved look much 
like those of other compressible convection simulations, 



dt 

dt 

djpE) 
dt 



It +v-(pv) = 0, (1) 



+V-(pvv) =-Vp + pg, (2) 

+V • [pvE] = -V • (pv) - K{T, p)VT 

+enuciXi,p, T) - €u{p, T) (3) 



where E is given by. 



E = e+-v\ (4) 

and K[T, p) is the radiative conductivity. Notice, however, that there are source and sink terms 
for nuclear energy release and neutrino losses. Thus, we must solve a set of continuity equations 
for composition as well, 

+ V • (pXiv) = p — - , (5) 



dt "^ ' ' '^K dt 



nuc 



where, except for the factor of p, the right side is a system of reaction network equations. This 
means that convection in oxygen shells involves energy production and loss which is directly in the 
flow. 



To our knowledge, the most physically similar system to have been examined previously with 
multi-dimensional simulations is the core helium flash ( |Cole &: Deupree 1980 , Deupree & Wallace 



1987, Deupree 1996). The oxygen shell differs in that (1) the reaction rates involved have even 
more extreme temperature dependencies (roughly T^^), (2) the equation of state is dominated 
by radiation and electron thermal pressures rather than degeneracy and thermal pressures, and 
(3) convection is driven by the competition between heating by nuclear burning and cooling by 
expansion and by neutrino emission (rather than cooling by radiative diffusion, as is the normal 
assumption in convective systems). 



3.1. Timescales 

The timescales of the problem portray this connection between convective mixing and 
energetics. Figure 1 shows the various timescales estimated from the MLT solution from the initial 
one-dimensional, hydrostatic stellar evolution model. The upper plot shows that the convective 
turnover time across a given radial grid zone, 

''cnv — ! V*-*/ 

is very much shorter than either the nuclear energy timescale, 

t = ^ (7) 

where Eth is the thermal component of the specific energy, or the Kelvin-Helmholtz timescale, 

iKH = J, (8) 

where fi if the gravitational potential energy and L the luminosity (radiative and convective). 
This means that local convective response is fast enough to control heating and cooling. However, 
using MLT velocities, the integrated convective turnover time across the whole oxygen convective 
region is long (10^ seconds, or 550 sound travel times), so that the time to move material across 
the whole region is approaching both the nuclear timescale at the base of the convective shell 
(10 seconds), and the neutrino Kelvin-Helmholtz timescale (10^^^ seconds) through the oxygen 
shell. This implies heterogeneity and significant nonspherical variation. The timescale for energy 
flow by radiative diffusion, 

™''" |v-(if(r,p)VT)| ^' 

is orders of magnitude larger than all of these timescales, and, thus, is not important in this 
stage of evolution. For this reason, we need not calculate the effects of radiative diffusion in these 
simulations. 



10 



The lower plot shows the various abundance timescales, for nucleus i, where 

dX, 



ti — Xi 



dt 



(10) 



The burning timescale of oxygen is of the same order as the fastest convective timescale across a 
given grid zone, but orders of magnitude less than the turnover time across the convective region. 
Notice that many other nuclei have short timescales as well. 

The magnitudes of these timescales lead us to several conclusions. First, composition will 
not be homogeneous. This will feed back into the structure via the energy source terms. Second, 
the system has strong fluctuations in space and in time. Because the thermal energy timescale 
is of the order of both the Kelvin-Helmholtz timescale and full convective mixing timescale, the 
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Fig. 1. — Timescales for Convection and Burning. Mixing and burn times are comparable. 
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dynamics of convection and evolution of the star will be strongly varying. This hardly satisfies the 
requirements of MLT which mandate that the dynamics of convection reach a statistical steady 
state. The energy supply changes as fast as convection can move energy across the region. 



3.2. Evolution of the Velocity Field 

Even at the beginning of our simulations, we see that MLT does not correctly represent 
the dynamics of convection in oxygen shells. Because the fluid transport timescales are initially 
longer than the oxygen burning time in the flame zone, the oxygen flame undergoes a initial 
thermal pulse, much like an AGB star (but faster), and a sound wave is emitted. This sound 
wave takes 18 seconds to cross the oxygen shell and effectively erases the MLT velocity field which 
we imposed initially. This allows the shell to establish dynamical structures from very low noise 
levels and reduces any worry that the initial perturbations somehow determine the final dynamical 
state. This pulse is not due to a faulty mapping of the initial model onto the 2D grid, but to 
an inconsistency between the assumed MLT physics of the convective shell, and the actual 2D 
hydrodynamics with heating and cooling. It is due to a thermal inconsistency involving assumed 
and actual convective cooling, not an error in overall hydrostatic balance. 

Figure 2 shows snapshots of velocity evolution after times of 100, 200, 300, and 400 seconds. 
Before about 100 seconds, the dynamics are dominated by rising and descending plumes, which are 
driven by simple buoyancy. Small eddies (scale length < pressure scale height) form at the base 
of the shell and form/merge like those in most two dimensional convection simulations ( porter & 



Woodward 1994]) . By 200 seconds (11 sound travel times across the original oxygen convective 



zone), the region is fully convective and buoyancy braking is readily apparent as a retarding force 
to upward motions, causing swirls. A large eddy dominates the upper two-thirds of the convective 
region, and one eddy dominates in the lower third. At 300 seconds, a single downward plume 
spans the upper 2/3 of the zone, and buoyancy still dominates the fiow dynamics. By 400 seconds, 
buoyancy braking affects the dynamics not only at the boundaries, but also throughout the zone. 

These simulations exhibit many of the generic traits of multidimensional compressible 
turbulent convection simulations. A comparison of these vortices with other two-dimensional 
convection simulations ( Hurlburt et al. 1986| , [Porter Sz Woodward 199^ ) shows that they are 



qualitatively similar in dimension (length/width ~ 1), and that a few large eddies dominate. 



Porter fc Woodward 1994 showed that this is to be expected, as two-dimensional simulations 
exhibit the merging of counterrotating eddies into larger ones. Earlier in our simulations, the 
lower 1/3 of the domain was populated by small eddies. 

Velocity magnitudes in these eddies are typically at Mach numbers of 0.1 to 0.2, which is 
more than a factor of ten higher than MLT would predict. The spatial extent of these eddies is on 
the order of 3 pressure scale heights. The simulation at this epoch is still time dependent so that 
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Fig. 2. — Snapshots of Velocity Fields, 100-400s. Large eddies dominate. 
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its structures should not be compared to the mixing length used in MLT. Only the radial velocity 
correlation length of a statistical steady state ( phan et al. 1982 , Chan fc Sofia 1989| , Porter & 



Woodward 1994| ) should be used as some measure of a mixing length. As we explain later, the 



dynamical coupling of the energy generation to the kinematical mixing, as well as the limited 
number of convective turnover times involved in the entire stage of shell oxygen burning, prevents 
a statistical steady state from ever being attained in this stellar evolution stage. 

In light of these results, we view the details of one dimensional stellar evolution calculations 
of the structure of the oxygen shell with skepticism. The smooth, steady flow assumed in those 
simulations appears to be a phenomenon which does not survive in more realistic calculations. 



3.3. Energy Fluxes 

The convective energy fluxes in this stage show both similarities and differences in comparison 
to previous compressible convection simulations. Figure 3 shows the azimuthally averaged 
enthalpy, kinetic, acoustic, and radiative fluxes after 400 seconds. The top panel shows the inner 
part of the O convective shell; the bottom panel gives the outer part, on different scales. Notice 
that the fluxes in the outer region are smaller in amplitude by more than a factor of 20. The 
enthalpy flux pictured here is negative through much of the Schwarzschild unstable region; MLT 
assumes it is positive. The formal convective boundaries are indicated as vertical dashed lines. 
The convective fluxes extend beyond these limits. 

As in previous convection simulations, we find a non-zero flux of energy to be carried in the 
form of sound waves, 

Fp = Vr ■ 6P, (11) 

and kinetic energy, 

i^fc = ^('pv • vvA (12) 

MLT conventions require that the net radial kinetic energy flux be zero and that pressure 
equilibrium is always in effect. However, compressible convection simulations (Hurlburt et alj 



1984| , |Chan et al. T982| , [Chan fc Sofia 198q , [Cattaneo et al. 1991| , Porter fc Woodward 199^ ), 



reveal a finite kinetic energy flux due to cool, dense material being compressed into plumes. Since 
these plumes occupy a filling fraction at constant radius which conserves mass flux relative to 
slower, upward moving material, the kinetic energy flux must be nonzero. In this particular model, 
the kinetic energy flux can reach as much as 20 % of the enthalpy flux, which is less than the 
peak values of 60 % seen in previous "box" simulations. One explanation for this difference might 
be that simulations of convection in a box are governed by the competition between convective 
instability and viscous dissipation by radiative (or thermal) diffusion. Here, convective motion 
competes with neutrino losses, as well as numerical dissipation. The estimated radiative flux is 
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at least 4 orders of magnitude less than the enthalpy flux at all points in the simulation. This 
fact was emphasized in early stellar evolution studies of massive stars (Arnett 1972a), where MLT 
theory yielded the result that the luminosity at the boundary of such convective zones would be 
zero. 

Figure 4 shows the enthalpy and kinetic fluxes and their sum at the top, middle, and bottom 
of the convective shell as a function of time. The dot-dashed curve is the kinetic flux, the dashed 
curve the convective flux, and the solid curve their sum. At no time during the simulation is the 
total luminosity equal to zero in the formally stable region. We note in passing that it is fortunate 
that neutrino losses govern the extent of this convection, since the numerical dissipation of the 
PPM method at this resolution ( Porter &: Woodward 199^ ), while yielding effective Reynolds' 
numbers ranging from 10^ at the base of the shell to 10^ at the top, also imply Prandtl numbers 
ranging from 200 at the base to 1000 at the top. The small numerical dissipation is still larger 
than the tiny radiative dissipation. This suggests that we would need at least a tenfold increase in 
resolution before the radiative diffusion could be quantitatively evaluated, if at this tiny level it 
were still important. 

An important component of the energy flow is the acoustic flux. This quantity measures the 
energy carried in the form of sound waves and is relevant to the kinetic energy flux in flows of 
significant Mach number ( [Hurlburt et al. 1984 ). In places the acoustic flux is a significant fraction 
of the kinetic flux. Thus, any attempt to model system,s such as these with im,plicit or anelastic 
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Fig. 3. — Pressure Contributions vs Radius. Thermal electron and radiation pressure dominate. 
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methods will miss the significant contribution to energy transport by sound waves. 

The enthalpy flux distribution is negative through much of the Schwarzschild unstable 
domain. This is because the neutrino losses are dominant globally, causing contraction. After the 
first sound crossing time, the simulation adjusts to a state representative of the physics involved. 
The negative enthalpy flux implies that the energy arising from the oxygen burning flame zone is 
insufficient to suppport the structure predicted by MLT, and energy is being drawn from above. 
Neutrino-cooled matter is contracting down on the flame zone. 

We find that there is a time dependance in many quantities throughout the duration of the 
simulation. For example, significant deviations in flux about an average value are seen in Figure 4. 
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Fig. 4. — Energy Fluxes at 400s. Acoustic (p), kinetic energy (k), enthalpy (c), and radiative (r) 
fluxes are shown. Enthalpy flux is often negative, and acoustic and kinetic fluxes are important. 
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However, there may be a trend toward to a roughly constant amphtude for each of the three 
regions after 200s. 

There are several reasons for the time variability in the simulation, and for it to continue until 
the time of core collapse. First of these are the similarities of both the nuclear and convective 
turnover timescales to the time left before core collapse. We estimate that there should be no 
more than 12 turnovers from the birth of the convective oxygen shell until core collapse, as 
estimated from our one-dimensional model. Previous box models of convection certainly required 
more turnovers before reaching a time- independent state ( Hurlburt et al. 1984 , plurlburt et al. 



19861) . The other reason for time dependence in the simulation has to do with the density contrast 
and the aspect ratio. Previous work stressed that for density contrasts of 21, a minimum aspect 
ratio (horizontal /vertical dimensions) of about 3 was necessary for a static state to be reached in 
80 turnover timescales ( [Hurlburt et al. 1984| ). In fact, since our simulation covers such a large 
spherical geometric area, the aspect ratio ranges from ~ 0.2 at the base to 7r/2 at the upper edge 
of the convective zone. Even for a 27r simulation in r and (p, rather than r and 9, the ratio is 
only four times this value. Simulations still should be performed to see if these trends still hold 
in arbitrary geometries (cylindrical -|- spherical) as the effects of adiabatic cooling on individual 
mass elements is more enhanced in a curvilinear coordinate system than a cartesian system. 
Nevertheless, assuming the trends from the cartesian simulations roughly hold true in a spherical 
system, and that a density contrast of about 200 exists across the convective zone, and the values 
of our aspect ratios, our simulations suggest that we can hardly expect this convective oxygen 
burning shell to ever reach a dynamical steady state before core collapse and that the application 
of mixing length theory in convective shell oxygen burning is unjustified. 



3.3.1. Penetration and Overshoot 

We now consider the presence of overshoot in these calculations. The structure below and 
above the unstable region is not changed much from the initial state, hence the designation of 
overshoot and not penetration. As Figure 3 showed, thermal energy is transported out of the stable 
region while kinetic energy is transported into the region, apparently not in quantities sufficient 
to alter the structure in these regions. However, the structure of the temperature and density at 
the lower boundary is significantly altered. Most of this has to do with the mixing timescale being 
slower than the local oxygen burning timescale, in which case the oxygen is depleted and Yg drops 
locally due to electron capture on the oxygen burning products. This, in turn causes changes in 
the entropy and composition, which then feeds back into the temperature and density. 

While these mass motions (Reynolds stresses) are not sufficient to alter the stellar structure 
significantly, they do allow mixing of composition from the convectively stable region into the 
unstable region. Figure 5 shows the composition profiles of ^^C. The ^^C abundance in this 
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oxygen-burning convective shell is purely the result of undershoot from the upper, ^^C-rich 
regions. A blob of ^^C rich matter sinks down through the oxygen shell, finally halted by its own 
rapid burning. This is a new effect, not seen in one-dimensional simulations. We will discuss the 
consequences below. 

Unlike the box simulations ( Purlburt et al. IQSC , Hurlburt et al. 1994), we find that overshoot 
is more vigorous at the upper rather than the lower edge of the unstable region. This seems to 
be related to the relative stability of our upper and lower zones ( [Hurlburt et al. 1994 ). The 
presence of a nuclear shell at the bottom of the zone leads to a value of 7 which changes faster as 
a function of position than it does at the top of the shell. This means that the relative stability 



Top 




50 100 150 200 250 300 350 400 





_ 








4x10" 


- 




' \ 


— 




- 


Middle 


'ij 1 


- 


2x10" 


- 




fl'lv\; 


- 





- 


,. — .A^ /• 


Ay • \ 1 ,• , 


^-^ 


: 




2x10" 


- 




'; ,'■' " ^ 


- 


4x10" 


'- 









200 250 300 350 400 




Fig. 5. — Energy Fluxes: time variability. Even with this low noise initial state, convection is well 
established after 200s, but highly variable. 
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of the bottom of the zone is higher than the top of the convective zone. 



3.3.2. Perturbations 



The early appearance of 7-hnes from decay of ^^Ni, as weh as its effects on the Hght curve 
of SN 1987A, left theorists with the idea that mixing occurred, so that density and/or velocity 
perturbations must exist in the progenitor models (Chan &: Lingenfelter 1987, Arnett 1988| , 
Kumagai et al. 1988, Pinto Sz Woosley 1988| ). As hydro dynamical simulations have shown ( [Arnett 
et al. 19891 , [Fryxell et al. 1991| , [Miiller et al. 1991| , [Herant fc Benz 1992| , Perant et al . 1992|, 



Hachisu et al. 199^ ), these perturbations facilitate the mixing of ^^Ni toward the surface through 



Rayleigh- Taylor instabilities. The exact details of the perturbations, such as scale distribution and 
amplitude, have not been self consistently examined, but linear analysis indicates that instabilities 
are likely at composition interfaces ( Benz &: Thielemann 1990| , Miiller et al. 199l| ). 



Because of the behavior of the equation of state, as seen in its linearized form: 
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the nonspherical perturbations are correlated, and show a strong temperature dependence due to 
the radiation pressure. 

Figure 6 shows the different contributions to the pressure, as a function of radius. Electron 
degeneracy and thermal ion pressure contribute little to the total. Roughly 60 % is from thermal 
electrons, while almost all of the remaining 40 % is radiation pressure. This simplifies the above 
expression to one of the form: 
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This expression is a good description of the fluctuations in the flows. 



(14) 



In Figure 7, we show the density perturbation structure in and around the oxygen shell at the 
end of the calculation. The perturbation contours range over ±10%. The largest perturbations 
are in the flame zone (the inner convective boundary), and at the outer convective boundary, as 
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may be seen in the right-hand panel. Notice that the thin flame zone, which lies near the leftmost 
arc in the left panel, has large perturbations in its small physical volume. The perturbations at 
the top of the convective zone have a much larger length scale, although their angular scales are 
similar to those in the flame zone. The "perturbation amplitudes" implied by the convection in 
the oxygen shell easily exceed the minimum needed to give significant mixing in supernova ejecta 
models (i.e. a few percent). Not surprisingly, the highest amplitude perturbations are found at 
the point where the linear analysis indicated: at the composition interfaces. 

This may be sufficient to explain the mixing of most of the ^^Ni observed in SN1987A, but 
apparently not that at highest velocity {v ~ 3,000 km s^^; 



see 



Arnett 199l| , jHerant fc Benz 19921) . 



The production of this high velocity material may be related to the passage of the explosion shock, 







Fig. 6.- 



12 



C Plume Snapshots. Such effects are not in ID calculations. 
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and explosive oxygen burning in this inhomogeneous shell ( [Arnett 1994 , Bazan and Arnett, in 
preparation). It might also be due to mixing outward by overshoot of ^^Ni just prior to explosion 
( Herant fc Benz 1992| ). This would require very high temperature oxygen burning; the silicon 
burning core has too large a neutron excess to produce ^^Ni. We have seen no evidence for such 
high temperature oxygen burning; probably the best chance would be at the very onset of core 
collapse. 

As with previous convection simulations ( Hurlburt et al. 1986 , Hurlburt et al. 1994| ), we 
conclude that gravity waves are responsible for the generation of the perturbations. A fourier 
analysis of the power in vertical velocities at different frequencies i^, was made for the time interval 
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Fig. 7. — Density Perturbations at 400s. Left: Contours to ±10%. Right: azimuthal averages 
(rms). 
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from 200 to 400s. Figure 8 shows the power as a function of frequency just above the convective 
zone (top panel), in the middle of the zone, and just below the convective zone (bottom panel); 
compare with Figure 4 of Hurlburt et al. 1986[ The solid, dashed, and dot-dashed lines were taken 



at 



113 

4' 2' 4 



of the angular domain. In the lower panel, we see the signature of a complex series of 
modes. The local Brunt- Vaisala frequency is Nbv = 1.55s~^, which is much higher than most 



of the power. As did [Hurlburt et al. 1986| , we interpret this as a time variation due to internal 
gravity waves, recalling that such waves have frequencies less than the buoyancy frequency. 
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Fig. 8. — Power Spectrum over 200-400s, for top, middle and bottom of the convective region. 
At the top, most of the power is below the Brunt- Vassla frequency, indicating the importance of 
acoustic and gravity modes. The spread in power over frequencies at the bottom (flame zone) is 
due to extreme time variability. 
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However, just above the convective zone (top panel), we find that the power is also at 
frequencies just below the local Brunt- Vaisala frequency, Nbv ~ 0.02s~^. The level of power here 



is larger, unlike Hurlburt et al. 198(: , where the power above the convection zone was smaller than 
the power just below the convection zone. The difference seems to be in the nature of the driving 
and damping. In our simulation, the driving energy flux from the oxygen burning flame zone 
varies by factors of two at the base of the Schwarzschild unstable region. In addition, the energy 
release rate in the ^^C-rich plume exceeds that from the oxygen flame. It has already been inferred 
that g-mode waves of high radial wave-number can appear at the interfaces between convectively 
stable and unstable regions ( [Press 1981 , Hurlburt et al. 1986| ), and in this case we seem to have 



preferential driving of these modes. 

While the appearance of g-mode wave interactions in this simulation is consistent other 



studies (Press 1981, Hurlburt et al. 1986| ), the details are different. Thermal diffusion is not 



present in this calculation, so that any damping of linear internal waves is the result of numerical 
dissipation (shocks). Pressure waves generated within the convective zone should have a constant 
horizontal wavenumber as their source. These internal waves are known to undergo internal 
reflection when N(r) = w, where N(r) is the Brunt- Vaissala frequency at radius r (Press 1981). 



We note that the detailed coupling of acoustic and gravity waves to the convection also depends 



on boundary conditions ( Hurlburt et al. 1994 ) 



The greatest power in the middle of the convective zone occurs at zero frequency (center 



panel), as found by Hurlburt et al. 1986 . However, despite our significantly better numerical 
resolution, our peaks are less well developed. We attribute this to the fact that our convection is 
much further from statistical relaxation. This is a physical feature of the epoch of oxygen shell 
burning, not a limitation of the simulation, as our previous discussion shows. 

Figure 9 shows the power spectrum of density perturbations versus angular separations A^, 
at 400s. Three radii are indicated: a solid line indicates a radius above the convective zone, a 
dot-dashed line the convective zone, and a dashed line a radius below the convective zone. There 
is a fall off at low angular separations, A0 < 0.1. The zoning is much finer than that, and will 
allow clumping on a much smaller scale {I/tlq ~ 0.008). This fall off could be taken as a measure 
of the scale length of the perturbations to eventually be compared with observations, and as such 
is an important quantity, but two dimensional simulations such as these are biased. Because the 
vortices are pinned on the grid — their angular momentum vectors stay perpendicular to the plane 
of computation — cascading to smaller scales is inhibited. Thus we expect clumping down at least 
to scales this size; this issue must be addressed with fully three dimensional simulations. 



3.4. Additional Convective Mixing 
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We mentioned that the asymmetry in the upflows and downflows has direct consequences 
for convective mixing. In Figure 6, we showed a time history of ^^C abundance. In the initial 
model, ^^C exists at the 10"^ abundance level in the oxygen convective zone, and up to the 0.4 
level above the zone. As rising material overshoots into the high ^^C abundance area, ^^C-rich 
matter is entrained. Eventually, cooling causes fast moving downward plumes to mix ^^C deep 
into the convective zone. These fingers terminate at the point where ^^C burns faster than it can 
be advected. Nuclear energy hot spots appear at these points, and by the end of the simulation 
the energy output from these hotspots are a factor of 30 higher than the oxygen flame zone itself. 
Even ^He, which exists still farther away from the MLT convection boundary, is brought down 
into the convective shell, albeit at a low level. 



Density Perturbation Power Spectrum: 460 x 128 jr/2 P, t = 400 s 




Fig. 9. — Power Spectrum of density perturbations versus angular separations A0, at 400s. 
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An interesting phenomenon is the mixing of higher Yg into the convective sheh. Figure 10 
shows the perturbation structure of the neutron excess r] = 1 — 2Y(. at the end of the calculation. 
The non-uniform distribution at both the upper and lower boundaries is readily apparent. Since 
the 'mass cut' between the proto-neutron star and supernova ejecta is expected to occur within 
this convective zone, the presence of a non-uniform distribution of Ye may alter core collapse 
( Burrows fc Lattimer 1985 , [Burrows &: Fryxell 1992 ). A larger effect could be from the coupling 
of the electron capture with hydrodynamics, which would change the entropy and the size of the 
core ( [Aufderheide 1993| , Arnett 1996 ). Furthermore, once the delineation between neutron star 
and ejecta is made, explosive nucleosynthesis in the ejecta will reflect the distribution of neutron 
excess rj. The differences of Arj ~ 0.002 that we find here will radically alter the details of the 
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Fig. 10. — Contours of neutron excess t]. Explosive nucleosynthesis yields are particularly sensitive 
to such variations. Twenty contours are plotted, linearly spaced from 7/ = 0to7/ = 2x 10"'^ . 
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quasi-statistical equilibrium (QSE) nuclear burning process ( [Truran et al. 1966 , Michaud & Fowler 



1972 , Woosley et al. 197^ ) so that a non- uniform distribution of isotopes is produced. This may 



be a welcome result, as one dimensional models of explosive burning in massive stars have usually 
produced too many neutron-rich isotopes (Hainebach et al. 1974, Hartmann et al. 1985| ), unless 



the mass cut was carefully chosen (e.g., Thielemann, Hashimoto, & Nomoto 1990). 



Most significantly, this extensive mixing by plumes may indicate that our concept of 
convection is too narrow to properly describe the neutrino cooled stages of massive star evolution. 
If this stage is so drastically altered, why not earlier stages as well? 



3.5. Stellar Structure and Evolution 

We have found that the results of our two-dimensional computations are inconsistent with ID 
hydrostatic MLT results in many ways. As yet the evolution has proceeded about 10 % of the 
time from oxygen shell ignition to core collapse (for the ID evolution of this initial model). What 
can we infer about 2D effects on evolution? 

In ID models, under MLT conventions, the 0-burning region is wholly contained within the 
convective zone. This determines the evolution: the convective region burns until it changes its 
fuel to ashes, and there is little ingestion of new fuel or loss of ashes. Thus, at the exhaustion 
of fuel in the convective region, this layer of ashes is added to the underlying core. Contraction 
follows, and a new shell generally ignites just above the old, now defunct convective shell, and the 



process is repeated. This behavior has long been known ( Rakavy, Shaviv, fc Zinamon 1967 , Arnett 



1969); see Arnett 1996 , §10.1 for details. Arnett 1972a noted its possibly non-physical nature 



contrasted it to the case of a radiatively cooled shell, and also considered an alternative in 
which the convective region both lost ashes and ingested new fuel, so that it burned out in mass 
coordinate as a "wave." However, we are aware of no published evolutionary sequences of neutrino 
cooled shells except those calculated with the implicit assumption that the first picture is correct. 
This distinction involves the nature of the boundaries of the convective region. 

In these 2D simulations, within the first 25 seconds of evolution time, the flame zone 'flashed' 
and changed its thermodynamic state from convective instability to stability, thereby removing 
itself from the convective zone, as MLT would define it. The resulting sound wave does not 
develop into a shock as it moves down the density gradient. No structure changes occur by the 
passage of this wave. 

The upper panel of Figure 11 shows the azimuthally averaged values of energy production, 
temperature, density, and convective speed, as a function of mass coordinate, for the initial and 
final state in this simulation (that is, at and 400s). The variables of the initial model are 
represented by solid lines; those of the final model by dashed and by dotted lines. The most 
obvious difference is in the convective speed. In the initial (MLT) model, the velocities are 
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significant only in the convective oxygen burning shell and the convective helium burning shell. 
They are separated by a stable region containing C, Ne, Mg, and mostly O. In the final model, the 
convective speed is larger, and does not go to zero in this intermediate region. While in MLT the 
convective speed is a slowly varying quantity, in 2D simulations it fluctuates strongly. 

This additional mixing has affected the composition structure, which is shown in the lower 
panel of Figure 11. This is especially noticeable in the range 2.5 < M^/Mq < 3.2. It is 
interesting that the He burning shell, which is broadened by radiative diffusion, widens more due 
to convection. 




Fig. 11. — Top: Azimuthally averaged e, p, and T, at and 400s. Bottom: same for composition. 
While the gross structure is preserved, the flame zone changes. The velocity field is qualitatively 
different from ID, allowing mixing across stable regions. 
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A second important difference is in the thin flame zone, which is strongly compressed in the 
mass coordinate, and therefore less obvious in Figure 11. Consider the temperature and density 
structure. Just above the inner boundary, the temperature has a small hole and the density a 
corresponding peak. At first sight, the run of both variables is almost idential between the initial 
and final models, but in fact the boundary has changed. The steep abundance gradient below 
the oxygen flame zone is becoming shallower, and the oxygen flame is no longer in the most 
convectively active region. Similarly, the outer edge of the temperature "hole" is smoothed out. 



orii4rpj: Heating rate ( B.ODE+ll, T.OSE+W) 




Fig. 12. — Closeup of Flame Region: Contours of the total rate of nuclear heating. 



Figure 12 shows a closeup of the heating rate in the flame region. There is one velocity vector 
for each 25 zones. The 20 density contours are logarithmically spaced between e = 8.0 x 10^^ and 



7.0 X lO^^erg g ^ s ^. The oxygen burning shell is burning at e < 1.3 



]^Qi3gpg g 1 s 1^ and so is 



barely discernable as an inner ledge. The oxygen burning is not uniform, having several peaks. 
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Most of the energy release comes from carbon and neon bm'ning, as blobs having X(^'^C) of a 
few percent burn vigorously. These conditions (T ~ 2.3 x lO^K and p ~ 8 x lO^g/cc), occurring 
on a hydrodynamic time scale, are typical of "explosive nucleosynthesis" of carbon and neon. 

Notice that the "hot spots" correspond to downward moving matter. The width of the whole 
flame zone is about 5 times wider that of just X(^^0), and much more uneven. The velocities 
in the fastest burning zones (both oxygen and carbon) are relatively small, so the ashes are not 
immediately swept out into the most active convection. 

This situation is drastically different from the conventional picture suggested by ID models. 
It causes worry that the 2D evolution would have diverged from the ID at a time even earlier than 
that for our initial model. We are currently attempting to push this sequence to later times (core 
collapse) . 



3.6. Dimensionality 

Two dimensional simulations such as these clearly are more realistic than one dimensional 
simulations with ad hoc prescriptions for convection and mixing. Hydrodynamic flow is simulated 
self-consistently, and heterogeneity in space and time may be examined. Nevertheless, the 
assumed "rotational" symmetry (not presumed due to rotation but simply enforced by limitations 
in computational resources), is simply not valid. Simulation of a 3D sector, for example by 
extending the 460x128 by 64 zones in 0, is feasible on a massively parallel supercomputer. We are 
currently continuing this calculation with a PVM version of PROMETHEUS, and will to extend 
the simulations to three dimensions by the time of publication. 

It is unclear at this time what physical state will be described by a three dimensional 
simulation. It is well known that the energy cascade in two dimensional turbulence (lower to 
higher scales) is opposite that in three dimensions (higher to lower scales). This effect can be seen 
quite easily in our and other PPM-based simulations ( porter fc Woodward 1994 ), where small 



eddies disappear in favor of larger scale eddies. A recent comparison of two- and three-dimensional 
simulations of convection in a proto- neutron star (Miiller &: Janka 1997) also show an increase 



of small scale features in the three-dimensional simulation relative to the two-dimensional 
simulations, with an overall lower kinetic energy in the three-dimensional case. Thus, the overall 
convective land composition fields are expected to be more homogeneous in a three-dimensional 
case. 

However, the effects of a third dimension on overshoot are not at all clear. This is critical for 
the oxygen shell evolution in light of what we mentioned previously concerning carbon brought 
into the oxygen shell by overshoot and the inhomogeneous energy release. On the one hand, 
we know that overshoot is the result of individual plumes creating gravity waves at convective 
stability/instability boundaries. Such plumes are known to be more energetic in three dimensions 
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than two due to a reduced drag effect (two dimensional 'plumes' are tori in realty) if starting from 
the same initial conditions. We know, though, that the initial underdensity which characterizes a 
plume is a function of the kinetic environment, which, from the supernova simulation mentioned 
above and turbulence studies, should be reduced from its two dimensonal analog. Just which effect 
will predominate requires a three dimensional simulation. More insight can be obtained from laser 
experiments for related hydrodynamic problems, where Rayleigh- Taylor instabilities in 2D and 3D 
are examined experimentally within the context of inertial confinement fusion ( [Remington et aL 



1995 ). Testing of PROMETHEUS against ICF codes and NOVA laser experiments is beginning 



QKane et al. 1996| ). 



4. CONCLUSIONS 

Direct calculation of the hydrodynamic behavior of oxygen burning shells indicates that the 



process is more interesting than previously supposed. Previous 2D work ( Arnett 1994 , Bazan fc 
Arnett 1994| ) is confirmed and extended. 



Composition is not homogeneous in the convective region. The system has strong fluctuations 
in both space and time. Mach numbers in the flow are typically 0.1 to 0.2; pressure fluctuations 
and acoustic fluxes are not negligible. Any attempt to model systems such as these with implicit 
or anelastic methods will miss the significant contribution of energy transport by sound waves. 
This convective flow will not reach a statistical steady state prior to core collapse. While the 
convective motions are not sufficient to alter the stellar structure significantly through their 
Reynolds stresses, they do allow mixing of composition from the convectively stable region into 
the unstable ones. The evolutionary implications of this have not yet been explored. 

In light of these results, we view with skepticism the details of one dimensional stellar 



evolution calculations for this stage (e.g., Woosley fc Weaver 1995 , Thielemann, Hashimoto, & 
Nomoto 1990|) . 



Prior to the core collapse, perturbations in density develop in the oxygen shell which are 
sufficiently large to "seed" hydrodynamic instabilities which will mix the "onion skin" composition 
of the presupernova. This occurs in precisely the region in which ^^Ni is explosively produced by 
oxygen burning behind the explosion shock. 

Although the active heating and cooling by nuclear burning and neutrino emission is novel, 
we find good agreement with previous multi-dimensional simulations of convection. Our results 
seem insensitive to numerical details. 

Enlightening discussions with Prof. Phillip Pinto and Willy Benz are gratefully acknowledged. 
Information from Nick Brummell about work in progress is sincerely appreciated. This work was 
supported in part by NASA grant NAGW-2450, and by NSF grant ASTRO-9015976. 
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A. NUMERICAL EXPERIMENTS 
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Fig. 13. — Velocity vectors for Models D (upper left), B (upper right), E (lower left), and C (lower 
right). See Table 1. 



Figure 13 shows the velocity vectors for several numerical experiments: Models D, B, E, and 
C from Table 1. The two top panels show a change in resolution by a factor of 2 in each linear 
dimension. The higher resolution case is to the left. In the high resolution model, velocity vectors 
are the average vector over 4 zones, while the vectors in the lower resolution models of 90 and 
120 degree extents are 2 and 3 zone averages, respectively. The flow patterns are qualitatively the 
same, and show some quantitative similarity as well, which leads us to believe that this resolution 
is the minimum needed for a 'converged' solution. The lower left panel shows the effect of widening 
the wedge from 90 to 120 degrees; we now have three vortices in the midlayer instead of two and 
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a half. A larger difference may be seen in the lower right panel, which used reflecting boundary 
conditions. The vortices are all contained within 90 degrees. Despite this, the qualitative features 
are still similar. 
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Fig. 14. — Abundance contours for C for models D, B, E, and C, as in Figure 13. See Table 1. 



Figure 14 shows the corresponding panels for contours of ^^C abundance. In all cases we have 
vigorous burning of an entrained plume of ^^C-rich matter. 
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